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Abstract 

The long-term dynamics of many dynamical systems evolve on an attracting, in- 
variant "slow manifold" that can be parameterized by a few observable variables. Yet 
a simulation using the full model of the problem requires initial values for all variables. 
Given a set of values for the observables parameterizing the slow manifold, one needs 
a procedure for finding the additional values such that the state is close to the slow 
manifold to some desired accuracy. We consider problems whose solution has a singu- 
lar perturbation expansion, although we do not know what it is nor have any way to 
compute it. We show in this paper that, under some conditions, computing the values 
of the remaining variables so that their (m + l)st time derivatives are zero provides 
an estimate of the unknown variables that is an mth-order approximation to a point 
on the slow manifold in sense to be defined. We then show how this criterion can be 
applied approximately when the system is defined by a legacy code rather than directly 
through closed form equations. 
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1 Introduction 



The derivation of reduced dynamic models for many chemical and physical pro- 
cesses hinges on the existence of a low-dimensional, attracting, invariant "slow mani- 
fold" characterizing the long-term process dynamics. This manifold is parameterized 
by a small number of system variables (or "observables" , functions of the system vari- 
ables): when the dynamics have approached this manifold, knowing these observables 
suffices to approximate the full system state. A reduced dynamic model for the evo- 
lution of the observables can then in principle be deduced; the resulting simplification 
in complexity and in size can be vital in understanding and modeling the full system 
behavior. This reduction has been the subject of intense study from the theoretical, 
practical and computational points of view. Low-dimensional center-unstable mani- 
folds are crucial in the study of normal forms and bifurcations in dynamical systems 
(e.g. [13]); the theory of Inertial Manifolds and Approximate Inertial Manifolds [6,31] 
has guided model reduction in dissipative Partial Differential Equations; the study 
of fast /slow dynamics in systems of Ordinary Differential Equations is the subject of 
geometric singular perturbation theory (e.g. [7]). On the modeling side, the Boden- 
stein "quasi-steady state" approximation has long been the basis for the reduction of 
complex chemical mechanisms, described by large sets of ODEs ( [2], see also the discus- 
sion by Turanyi [33]). More recently, an array of computational approaches have been 
proposed that bridge singular perturbation theory with large scale scientific computa- 
tion for such problems; they include the Computational Singular Perturbation (CSP) 
approach of Lam and Goussis [21-23], and the Intrinsic Low-Dimensional Manifold 
approach of Maas and Pope [24]. The mathematical underpinnings of these methods 
have also been studied ( [15, 34]). Lower-dimensional manifolds arise naturally also 
in the context of differential-algebraic equations, where the initialization problem has 
attracted considerable attention (e.g. [3]). 

Remarkably, the same concept of separation of time scales and low-dimensional 
long-term dynamics underpins the derivation of "effectively simple" descriptions of 
complex systems. In this context, the detailed model is a collection of agents (molecules, 
cells, individuals) interacting with each other and their environment; the entire distri- 
bution of these agents that evolves through atomistic or stochastic dynamic rules. In 
many problems of practical interest it is possible to write macroscopic equations for 
the large-scale, coarse-grained dynamics of the evolving distribution in terms of only 
a few quantities, typically lower moments of the evolving distribution. In the case of 
isothermal Newtonian flow, for example, we can write closed evolution equations for the 
density and the momentum fields, the zeroth and first moments of the particle distribu- 
tion over velocities. This is again a singularly perturbed problem; only in this case the 
higher moments of the evolving distribution have become quickly slaved to the lower 
ones (in this case stresses, after a few collisions, have become functionals of velocity 
gradients). Newton's law of viscosity therefore represents a similar type of "slow mani- 
fold" as in the ODE case discussed above - fast variables (stresses) become functionals 
of velocity gradients, and the slow manifold is embodied in the closure: Newton's law 
of viscosity. The use of slow manifolds in non-equilibrium thermodynamics, and more 
generally in the study of complex systems, is also a subject of intense current research 
(see the work of Gorban, Karlin, Oettinger and coworkers [11,12], as well as [18,19]) 
for some recent computational studies. In this context, one has an "inner simulator" at 
the microscopic or stochastic level for evolving the detailed distributions; a separation 
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of time scales does not arise at this level, but rather at the level of the evolution of the 
statistics or moments of these distributions. Typically the lower moments of the evolv- 
ing distributions parameterize the slow manifold, while the higher moments quickly 
evolve to functionals of the lower ones. Since we do not have explicit formulas for 
the equations at the coarse-grained, macroscopic level, the following interesting ques- 
tion arises: can we benefit from singular perturbation, when no closed form evolution 
equations are available, and the only tool at our disposal is a "black box" dynamic 
simulator of the detailed problem ? This is the problem we will address in this paper. 
We will assume that we are given an evolutionary system which can be described 

by 

u = p(u,v) (1) 
v' = q(u,v) (2) 

where prime designates differentiation w.r.t. t, the dimensions of u and v are N u and N v 
respectively, and values, n(0), are specified only for u. In the case of a legacy dynamic 
code, we may not even be given the formulas for these equations explicitly; instead, 
we may be given a time-stepper of the above system as a black-box subroutine: a code 
that, provided an initial condition for u and v, will return an accurate approximation of 
u and v after a time interval (a reporting horizon) AT. We wish to find v(0) so that the 
solution is "close to a slow manifold." This statement is deliberately vague because in 
practice we are proceeding on the belief that there exists a slow, attracting, invariant 
manifold that can be parameterized by u and that the variables v, in some sense, 
"contain" the fast variables so that their values on the slow manifold can be computed 
from the values of u at any time. (Note that we do not need to know which are the slow 
variables, only to be able to identify a set of variables sufficient to parameterize the slow 
manifold.) This implies that the manifold is the graph of a function v = v(u) over the 
observables u. As we proceed, we will make these statements more precise. However, 
we will not cast them into the form of theorems because, even though this is possible, 
the conditions for the application of the theorems would be essentially impossible to 
verify in all but trivial problems. As with many numerical methods, the primary test 
of applicability is in the actual application. 

We assume that the system can be expressed in terms of other variables, x and y, 
of the same dimension as u and v respectively, and that in terms of x and y the system 
can be written in the usual singular perturbation form: 

x' = f(x,y) (3) 
ey' = g(x,y) (4) 

where x and y are also of dimension N u and N v , respectively. We will assume that their 
initial values are specified as x(0) and 2/(0), independent of e. The singular perturbation 
parameter, e, is associated with the gap, or ratio, between the "active" (slow) eigenval- 
ues and the rightmost of the negative, inactive eigenvalues of the linearized problem, 
locally. We stress that we do not assume that we know how to express the equations 
in the form of eqs © and (J3J) nor do we have any knowledge of the transformation 

u = u(x,y), x(0, e) = x(0) (5) 
v = v(x,y), 2/(0, e) = y(0) (6) 
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other than that we assume that it is well-conditioned and does not have large deriva- 
tives. We will consider singular perturbation expansions in e, even though the param- 
eter is not identified (and cannot be varied). The functions / and g could also involve 
the parameter e, but that serves little in this presentation other than to complicate the 
algebra. 

The standard singular perturbation expansion for the solution of eqs © and (jlj) 
takes the form: 

oo oo 

x(t,e) = ^aW + ^eU) (7) 

n=0 n=0 

oo oo 

y(t,e) = £ e n Y n (t) + £ e n Vn (r) (8) 

n=0 n=0 

This involves an outer solution (X(t),Y(t)), that is smooth in the sense that its time 
derivatives are modest, and an inner solution, (£(t), tj(t)), that captures the fast bound- 
ary layer where the solution typically changes like e~ */ e = e~ T . Both outer and inner 
solutions are expressed as power series in e, the latter as a function of the stretched 
time t. The inner solution is fast in the sense that each differentiation by t introduces 
a multiplier of 1/e. 

We define the slow manifold as the manifold that contains all solutions of eqs J3J) 
and Q of the form eqs (J7J) and (JHJ) with the inner solution identically zero. This is an 
invariant manifold. We say that a solution is an mth-order approximation to a slow 
manifold solution if it has the form in eqs (jJJ) and (JHJ with the first m + 1 terms of the 
inner solution expansion identically zero. A point (u, v) is an m-th order approximation 
to the slow manifold, or m-th order close to the slow manifold, if it lies on an m-th 
order approximate solution. 

We want to stress that we are not proposing a technique for finding the singular 
perturbation expansion. Rather, we are using the ideas as a scaffold for the theoretical 
justification of the proposed computational method. It is possible that the method 
will provide answers even when the singular perturbation expansions do not converge 
(although in that case we have no justification other than an intuitive one). The 
procedure we propose for finding the v values that are close to the slow manifold 
given the u values is to find values of v that approximately solves the N v dimensional 
non-linear equation 

= (9) 

which we call "the [(m + l)-st] derivative condition." (Compare this with the "bounded 
derivative principle" [20] which requires the first m time derivatives to be of order 1. 
That condition can be applied to problems with fast oscillating solutions. Ours can 
not, but it is simpler to apply to the types of problems we consider.) 

Note that eq. (JUJ) is a local condition - that is, it is applied at a single time - which we 
will take to be t = to = - to determine a value of v corresponding to a given value of 
u. A solution of eqs Q and © starting from these values of u and v will not satisfy eq. 
© for t > - but we do expect the solution to be close to the slow manifold. Intuitively 
the condition in eq. @ finds a point close to the slow manifold because differentiation 
"amplifies" rapidly varying components more than slowly varying components, so eq. 
(JHJ) seeks a region where the fast components are small. We will suggest ways in which 
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eq. © can be solved approximately in practical codes, even those based on a legacy 
code for the integration of eqs (JTJ) and @. 

While the approach will be presented and implemented for singularly perturbed sets 
of ODEs, the "legacy code" formulation is appropriate also in cases where the inner 
simulator is not a differential equation solver, but rather a microscopic / stochastic sim- 
ulator. In the equation free approach we have been developing for the computer-assisted 
study of certain classes of complex systems, the variables u correspond to macroscopic 
observables of a microscopic simulation (typically, moments of a stochastically or de- 
terministically evolving distribution). In this case v corresponds to statistics of the 
evolving distribution (e.g. higher moments) that become quickly slaved to (become 
functionals of) the observables u; thus the analogy of a slow manifold persists in mo- 
ments space for the evolving distribution. 

The paper is organized as follows: In the next Section we will outline the theoretical 
basis for the derivative condition. In Section 3 we discuss ways in which the derivative 
condition can be approximately applied as a difference condition and used with legacy 
codes. Section 4 presents some simple examples of its application. We conclude with 
a brief summary and outline of the scope of the method and some of the challenges we 
expect to arise in its wider application. 

2 Theoretical Basis 

In this section we will show that the application of the condition in eq. © will 
lead to a mth-order approximation to the slow manifold under suitable smoothness 
and smallness conditions. We will start by sketching the parts of singular perturbation 
expansion theory that we need by paraphrasing the presentation in O'Malley's mono- 
graph [26], particularly pp 46-52. We will use the same notation to make it easier for 
the reader who wishes to get more detail from that book. 

In the following we will write X n , Y n , £ n , and r) n to mean X n (t), Y n (t), £, n (T~), 
and %(r). We recall from [26] that the t and r dependencies are treated separately, 
and that the terms in the outer expansion, {X n ,Y n }, are obtained term by term by 
substituting eqs and (jHJ) into eqs ® and and equating each outer term in e n 
to zero, starting with n = 0. For n = we obtain the DAE: 

X f = f(X ,Y ), X (0)=x(0) (10) 
= g(X ,Y ) (11) 

The existence of a smooth solution of this equation for any x(0) requires the assumption 
that g y is non-singular. (The existence of asymptotic expansions for the inner and outer 
components requires the stronger assumption that g y is a stable matrix.) Then, Yq is 
specified uniquely in terms of Xq by eq. say as 

Yq = <t>(X Q ). 

The nth term in the power series yields the DAE 

X' n = f x (X ,Y Q )X n + f y (X ,Y )Y n + f n (12) 
= 9x (X ,Y )X n + g y (X ,Y )Y n + g n (13) 
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where the f n and g n are defined in terms of earlier terms in the outer expansion, 
{Xj,Yj},j = 0, • • • ,n — 1. The initial condition, Xj(0), has to be specified. Since we 
have set X(0) = x(0), Xj(0) is obtained from eq. (JTJ) by requiring that the e J term 
vanishes at t = 0, which gives 

Xj(0) = -£-1(0). (14) 

We are most interested in the way in which the inner terms are defined, since we 
wish to annihilate the first m + 1 of these to get an mth-order approximation. These 
are obtained by considering the change in the inner terms as r varies for arbitrarily 
small e, in other words, with t = and the outer solutions fixed at their initial values. 
Following [26] we consider terms in successive powers of e and find that the e° terms 
satisfy: 

£ = /(x(o),y (o) + %)-/(x(o),y (o)) (15) 

% = g(x(0),Y (0)+Tto)-g(x(0),Y (p)) (16) 

where a dot represents differentiation w.r.t. r = t/e. If we have an initial value for 
7?o(0) we can solve eq. (fT6|) for 770. Eq. eq-inOl gives £0 as an indefinite integral so it 
is determined by specifying £0 at any point. This is normally done at r = 00, in other 
words, at the end of the boundary layer. However, in our development here we will 
be showing that rjj and £j are identically zero for j < m, so we will actually choose 
0(0) = (so that we also have X, + i(0) = from eq. <fT3j) ). Subsequent inner terms 
satisfy 

L = fy(x(0),Y O {0)+ VO )Vn + fn (17) 

7?n = g y {x(0),Y (0) +T] )r] n + g n (18) 

where f n and g n are functions of the earlier terms, {£7,77,}, j = 0, ■■■,n — 1. In 
particular, if all of these terms are zero, then f n and g n are zero. Eq. eq-inn2 can 
be solved if an initial value is known for rj n (Q). Once again, eq. Q17|) gives £ n as an 
indefinite integral. 

In the following we are going to show, one by one, that rjj(0) = for j = 0, 1, • • • ,m, 
and that we can then choose £j(0) = 0. Note that once we have shown that 7/0 (0) = 
then eqs (|15|) and (|16[) indicate that 770(7") = and that £o( T ) is constant, which we can 
make zero by choosing £o(0) = 0. Then it follows that /1 and g\ are identically zero. 
Repeating this argument, we see that if 7?j(0) = 0,^ (0) = 0, j = 0, ■ • • ,m then £j(t) 
and ?7j (t) are identically zero for j < m. This provides a solution that is an mth-order 
approximation . 

Now we return to the original problem phrased in terms of u and v. If we knew 
the transformation to x and y we would do better to work in that space, but our 
assumption is that, although a transformation exists, it is unknown and we have to 
work with u and v. We want to show that if the (m + l)st derivative of v is zero, then 
the point is mth-order close to the slow manifold. All terms below are evaluated at 
t = 0orr = 0- the time at which we are attempting to solve eq. © • We will simplify 
the notation and write rjj for 77^ (0) and similarly for other terms in the following. We 
have from eq. ((HI 

d m+l v d m+l y 

— — rr- = v v — — - + other terms 19 

dt m+1 y dt m+1 
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where the other terms involve either partial derivatives of v w.r.t. x and/or multiple 
derivatives of v w.r.t. y and products of derivatives of y. 

Substituting from eq. jSJ) into eq. (|19j) and extracting the lowest order term in e 
(e-™- 1 ) we find that at t = 

r_ " = e~ m -\ , '*? + other terms + 0(e- ,n ) (20) 

where now the other terms include products of a higher-order partial derivative of 
v w.r.t. y with more than one r-derivative of ijo such that the sum of the levels of 
differentiation is m + 1, that is, terms like 



d k 7] d m+1 ~ k T] 
v vy dr k dT m+i-k 



d m+1 v 
dt m+1 



9y 1 Vo + other terms (21) 



+ 0(e" m ) (22) 



and terms with higher partial derivatives and more derivatives of rjQ in the product. 
Note that whenever r/o appears, it is always differentiated w.r.t. r at least once. Also 
note that we do not get any terms involving £o because of the additional e appearing 
in front of the inner solution expansion for x in eq. 

Now we use eq. ()16(l to find the higher-order derivatives of r/o w.r.t. r. We get for 
p > 1 

d p Vo = nP _i i 
drP 

where the other terms involve ?)q with j > 1. Substituting eq. (|21j) in eq. (|20() we 
arrive at 

m+l 

v y g™Vo + v ^ z rf 
i= 2 

where the notation v z g z stands for sums of products of various partial derivatives of v 
and g. Equating the leading term of the right-hand side of eq. (|22j) to zero, we now 
have a polynomial equation for t)o as 

m+l 

VyQym + J2 v z9*4 = (23) 

One solution of this is 

?)o = 

and it is an isolated root as long as v y g y is non-singular. Since we have assumed that 
g y is a stable matrix for the existence of a singularly perturbed solution (and hence a 
slow manifold) and that v in some sense spans the fast variables (meaning that v y is 
non singular), this is no problem. If all other partial derivatives involved in eq. (j2.3[) 
are "of order one" then other solutions are also of order one - i.e., well separated from 
the zero solution. We will delay discussion of how to avoid the "wrong" solutions for 
the moment, and assume that we find the zero solution. (If the problem is linear, 
these other terms are null, so there are no other solutions, and it is only in the case of 
high non-linearity when the partial derivatives are large that these other solutions can 
become small and cause problems.) 

If 7)0 = then eq. (|TH|) tells us that rjo = because we have assumed that g y is a 
stable matrix (in the domain of interest). This immediately implies that £0 = (or is 
a constant that can be absorbed into the outer solution, thus making £0 = 0). 
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As discussed following eq. I|18|). the vanishing of rj and £o means that the last 
terms of eqs (|17f) and (|18|) are zero for n = 1, making them look similar to eqs (|15|) 
and ()16() . Therefore the above argument can now be applied to show that r/i and £1 
are zero. This argument can be repeated for higher-order terms as long as the power 
of e in eq. (|22|) remains negative, in other words, until we have shown that 

= €j = °> i = 0, • • • , m 

Note that when we have made the (m + l)st derivative zero, the lower order deriva- 
tives will not be zero, or even small. This is because a small movement away from the 
slow manifold can make large changes in the derivatives of the inner solution. However, 
the difference between successive v values as we make m successively larger is small - 
of order e m+ as we go from the m-th to the m + 1-st derivative condition as the v 
values are converging to the slow manifold. Hence one way to solve for zeros of high- 
order derivatives would be to start by finding the zeros of the first order derivative, 
then repeating for successively higher-order derivatives, each step requiring smaller 
and smaller changes to v, until we have found the zeros of the (m + l)st derivatives 
using whatever computational process is appropriate. (The computational process is 
addressed in the next section.) This procedure helps address the issue of finding the 
smallest of multiple roots of eq. (|23|) since, for m = 0, there is only one root so that 
the iteration for m = 1 and larger m starts with a good approximation. If the zero 
root is well separated from the others, we will converge to it. 



3 Practical Implementation 

It is often not practical to work with higher-order derivatives of a differential equa- 
tion, either because they are algebraically complicated or because the equations are 
defined by a "legacy code" - that is, as an implementation of a step-by-step integrator 
that effectively cannot be change or analyzed. (The same would be true if part of the 
derivative calculations involved table look-up functions that could be difficult to differ- 
entiate.) Therefore, we are interested in methods that do not require direct access to 
the mathematical functions constituting the differential equation. The same rationale 
applies when the "inner simulator" simulates the system at a different level (i.e. in the 
form of an evolving microscopic or stochastic distribution). In this case we only have a 
time-T map for the macroscopic observables, that results from running the microscopic 
simulator and monitoring the evolution of the observables (e.g. particle densities) in 
time [10,32]. 

The obvious alternative is to use a forward difference approximation to the deriva- 
tive, noting that if 

A rn+1 v{t) = A m v(t + H)- A m v(t) (24) 



with A°v(t) =v(t), then 

A m+1 v(t) = H m+1 ^——^- + 0(H m+2 ) 

v i dt m+i y > 

It turns out that there is a straightforward way to implement a functional iteration 
to find a zero of the {m + l)-st forward difference, even if we only have a " black box" 
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code that integrates eqs Q and (J2J). Suppose we have operators, <f> and 6, that, given 
values of u(t) and v(t), yield approximations to + H) and u(i + H), namely 

u(t + H) ^ (f>(u(t),v(t)) (25) 

u(t + H)»e(u(t),v(t)) (26) 

Letting t n = to + nil and applying eqs (|25j) and (|2f)|) m + 1 times starting from t = to 
we can compute approximations 

Uj+x = (f>(uj,Vj) ps u(ij+i) (27) 

Vj+i = 0(uj,Vj) ps (28) 

for j = 0, • • • , m. The functional iteration to find a t>o for a given iio such that the 
{m + l)-st difference is zero consists of the following steps: 

1. Start with a given uo and a guess of vq. 

2. Set the iteration number p = 1. 

3. Set the current iterate = vq 

4. Apply eqs ([27jl and p8|) m+1 times starting from uq, to generate , , • • • , i 

5. Compute 5 = (-l) m A m+1 v { p) 

6. If 5 is small, the iteration has converged. 

7. If 5 is not small, 

v W= v M + 6 (29) 

8. Increment p and return to step 4. 

If the black box integrator provides a good integration (that is, it does not introduce 
spurious oscillations due to near instability) and H is small enough, this process will 
converge on a zero of the difference. 

As an illustration we consider the case m = and assume that the integrator is 
simply forward Euler with a step size such that its product with the magnitude of the 
largest eigenvalue is less than one. We see that the process for v consists of computing 

dv 

5 = Vl - Vq « H—(to). 

at 

If this is insufficiently small, we replace vo with vq + 5 = v\. This is just the stationary 
projection process used in [8] and is related to the "reverse time" projective integration 
method in [9] . In the case of m = 1 we compute 5 as 

5 = —V2 + 2vi - Vq 

If this is insufficiently small, we add it to Vq to get a new vq given by 

Vq = 2vi - V 2 

This is precisely the linear interpolant through v\ and V2 back to the starting point. 
It is illustrated in Figure ^ (This Figure may be a little confusing because it is 
drawn in the u — v plane to emphasize that u is being held constant from iteration 
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to iteration. The backwards interpolant, however, is really taking place in the t — v 
plane.) The general iteration is equivalent to the obvious extension of that: an rath 
degree interpolant is passed through v%, V2, ■ ■ ■ , v m+ i to compute a new approximation 
to vq. This is conveniently done using differences as described above. 

Why does this iteration eq. ()29j) converge for small H under reasonable circum- 
stances? Intuitively, the forward integration exponentially damps the fast components 
more than they are amplified by the polynomial extrapolation backwards. In more 
mathematical terms, the iteration takes the form 

v% ew = v + 6 (30) 

If 85/dvQ is negative definite and small, this will converge. We have 

d5/dv = -f-m m+1 ^?^- + 0(H m+2 ) 

ov 

The term 

dv 

is dominated by 

m+l m+l 
fc Uy 

in powers of e. Since we have assumed that g y is a negative definite matrix for the 
existence of a singular perturbation expansion, convergence follows for small enough 
H. 

In the above discussion we have used the attractivity of the manifold and, in effect, 
successive substitution in order to converge to a fixed point of our mapping. One can 
accelerate this computation by using fixed point algorithms, like Newton's method; 
clearly, since no equations and Jacobians are available, the problem lends itself to 
matrix-free fixed point implementations like the Recursive Projection Method by Shroff 
and Keller [28] or Newton-Krylov implementations (see [16] and, for a GMRES-based 
implementation for timesteppers [17]). 



4 Examples 

We will consider three examples to illustrate the method. The Michaelis-Menten 
enzyme kinetics model is a classic example for singular perturbation (see [1] for a brief 
discussion on the history of the model and its analysis). Since it is a simple system 
we can make a direct comparison with an easily computable singular perturbation 
expansion. The second example is a realistic five-dimensional chemical reaction system 
with a one dimensional slow manifold. As in many real examples, we do not know the 
slow manifold, and can only show "plausibility" of our solution. The final example 
is a contrived five-dimensional non-linear system with a known two-dimensional slow 
manifold so that we can compute the "errors" as the distance from the slow manifold. 

4.1 Michaelis-Menten equation 

This is given in a singular perturbation form in [26] as 

x = —x + (x + k — \)y (31) 
ey' = x — {x + n)y (32) 
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Table 1: Michaelis-Menten example with e = 0.1, h = e/100 and n — 1. 



m 


Computed ?/ 


Asymptotic y 





0.500000000 


0.500000000 


1 


0.503049486 


0.503125000 


2 


0.503031986 


0.503027344 


3 


0.503031924 




4 


0.503031924 





with x(0) = 1. For some simple tests we will take k = 1, X = 0.5, and e = 0.1 and 
0.01. (These are larger figures than typical for reactions, but we wish to show that the 
method works even for problems with a relative small gap, and also to have problems 
where the higher-order initializations are visibly different from the lowest order ones.) 
The first few terms of the outer solution are 

t- kXx kXx(2kX — 3xX — xk — k 2 ) 9 ^. q% , , 

+ 7Z i ,^4 g + 7ZrT^7 ~t + °( e )• (33) 



x + K (x + k) 4 (x + k) 

In the following tests, we implemented the operators and 9 in eqs (|25|) and (|26j) 
using n steps of forward Euler with step size h. In all cases eq. (|5U)) was iterated until 
6 was less than 10~ 14 (which is rather excessive, but we didn't want any errors from 
premature termination to color the results). We ran with m = 0, • • • , 4 and compared 
them with the first m + 1 terms of eq. ()33|) through m = 2. For e = 0.1, the results are 
shown in Tabled with h = e/100 and n = 1, and in Table [21 with h = e/10 and n = 4. 
The tables show the computed approximation to the slow manifold, y(0), for x(0) = 1. 
The tables also show the first m terms of the asymptotic expansion for m < 2. As can 
be seen, the discrepancies are larger when H = nh is larger. 

A larger "integration time horizon" H gives more rapid convergence. However, this 
means than the difference estimate is a less accurate approximation to the (m + 1)- 
st derivative. Our theory shows that making the (m + l)-st derivative zero puts the 
solution on an m-th order approximation to the slow manifold. Any error in the 
derivative estimate creates an error of the same order in the solution, so we should 
choose the time horizon so that the errors in the derivative estimate are of the same 
order as those we are willing to tolerate in the approximation to the slow manifold. 
This suggests that an initial approximation could be calculated with a larger H (and 
small m) to get faster convergence, and then H could be reduced as m is increased to 
increase the accuracy. 

The next case, shown in Table 01 uses an e smaller by a factor of 10. Because 
the higher-order terms in e are now smaller, the agreement with the terms in the 
asymptotic expansion is better. (However, we are really interested in the agreement 
of the computed vo with the v on the slow manifold, rather than with the ability to 
match the first few terms in the asymptotic expansion.) 

The cases above "constrained" the derivative of the singularly perturbed "fast" 
variable, y. Usually we cannot isolate this variable. To see the effect of having a 
different variable set, we transform x and y into 

u = x + y (34) 
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Figure 1: Iterative process for m = 1. 



Table 2: Michaelis-Menten example with e = 0.1, h = e/10 and n = 4. 

m Computed y Asymptotic y 

~0 0.498886090 0.500000000 

1 0.503067929 0.503125000 

2 0.503035446 0.503027344 

3 0.503035098 

4 0.503035128 



Table 3: Michaelis-Menten example with e = 0.01, h = e/100 and n—1. 



m 


Computed y 


Asymptotic y 





0.500000000 


0.500000000 


1 


0.500311725 


0.500312500 


2 


0.500311533 


0.500311523 
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Table 4: Michaelis-Menten example with e = 0.1, h = e/10 and n — 4. 
m x y ytrue 

~0 0.98825957 0.51174043 0.50011008 

1 0.99743598 0.50256402 0.50239316 

2 0.99756721 0.50243279 0.50242566 



Table 5: Michaelis-Menten example with e = 0.01, h = e/10 and n — 4. 



m 




2/ 


ytrue 





0.99874363 


0.50125637 


0.49999762 


1 


0.99974927 


0.50025073 


0.50024891 


2 


0.99975069 


0.50024931 


0.50024927 



v = y - x, (35) 

and work with u and v assuming that we do not know this transformation. (We chose 
this transformation because in some sense it puts equal parts of the slow and fast 
variables, x and y, in u and v, illustrating the fact that we need only know variables 
that parameterize the slow manifold (u in this case), not ones that are in some sense 
dominated by the slow manifold.) We assumed that we were given a value of u, 
u(0) = 1.5, and computed the approximation to v(0) using our method. This was run 
with h = e/10 and n = 4 for e = 0.1 in Table H and for e = 0.01 in Table El The 
tables show the corresponding x and y values derived from u = 1.5 and the computed 
v(0). The column labeled ytrue gives the first three terms of the asymptotic expansion 
of y given the x shown in the first column. Now we can see the significant error that 
the m = 1 case gives rise to when we "come at the slow manifold from a different 
angle." However, the higher-order approximations yield a good approximation to the 
slow manifold. 

4.2 A Simplified Hydrogen-Oxygen Reaction System 

We will use an example from Lam and Goussis [23]. It contains seven radicals, O2, 
H, OH, O, H2, H2O, and HO2 which we will group in that order as the vector y. The 
differential equations are: 



dm 

dt 

dm 

dt 

dm 

dt 



-hfyw2 + hby-m + ^sym - ^5/2/12/2 

-hjym + k u y3y4 + ^2/2/42/5 - £252/22/3 + £3/2/32/5 - %2/22/6 - ^£5/2/12/2 

^1/2/12/2 - &162/32/4 + £2/2/42/5 + £262/22/3 

-£3/2/32/5 + £3b2/22/6 - £4/2/32/7 - 2£ 8 /2/I + 2k 8 bym 



= £i/2/i2/2 - £lf>2/32/4 - £2/2/42/5 + £262/22/3 + £8/2/1 ~ £862/42/6 

-Jr = -£2/2/42/5 + £262/22/3 - £3/2/32/5 + £362/22/6 
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Table 6: Reaction Rates for Hydrogen Oxygen System 
% kif ktf, 



1 1.0136 x 10 12 1.1007 x 10 13 

2 3.5699 x 10 12 3.2105 x 10 12 

3 4.7430 x 10 12 1.8240 x 10 11 

4 6.0000 x 10 13 

5 6.2868 x 10 15 

8 6.5325 x 10 12 3.1906 x 10 11 



Ta ble 7: Eigenval ues 
-2.5 x 10 6 
-1.4 x 10 6 
-4.0 x 10 4 
-8.3 x 10 3 
-4.0 x 10~ 3 



—j^- = h f y 3 y 5 - hbV2VQ + h f y 3 y 7 + hfvi - hbVAye 
^ = -k if y 3 y 7 + /x/c 5 /yi?/2 

The values of the coefficients are taken from the cited paper and shown in Table H3 
The parameter \i represents pressure and is 4.5 x 10 -4 . The differential system has two 
constants of integration representing mass balance for oxygen and hydrogen atoms, so 
it is really a five-dimensional system. The eigenvalues of a local linearization in the 
region of operation for this example are approximately as shown in Table [7| From 
these we see that after an interval of order one millisecond, the system is effectively 
one dimensional. 

In the test below we chose 7/5 (O2) as the observable variable. The system was run 
from the starting conditions given in [23] until t = 6.41 x 10 -4 (one of the reporting 
times in their paper). Then our process was applied, fixing y§ to its current value, 
and choosing all other variables so that their (m + l)-st forward difference is zero. 
To emulate a "legacy code" situation, we integrated the equations using a standard 
integrator (ode23s in MATLAB) over m + 1 intervals of length H. In this example, 
we used H = 1 x 10~ 5 . The relative and absolute error tolerances for ode23s were set 
to 10~ 12 and 10 -14 , respectively. To maintain the mass balance relationship, radical 
concentrations of H and OH are computed directly from the mass balance relations. 
(These two were chosen because their concentrations remain reasonably non-zero. If a 
radical whose concentration gets close to zero is used, there is some danger of roundoff 
errors causing the concentration to become negative. This will often make the system 
unstable - as well as physically unrealistic.) 

The procedure was run with m = 0, 1, • • • , 4. The starting values of the concentra- 
tions were as shown in Table |H1 (These are given for the sake of completeness should 
anyone wish to compare with our results.) The constrained results for each m differ 
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Table 8: Radical Concentrations Prior to Constraining to Slow Manifold 



7/1 


4 2783465727 x 10" 


-13 


1)2 


3.9878034748 x 10" 


-8 


2/3 


1.3883748623 x 10" 


-10 


2/4 


1.1300067412 x 10" 


-11 


2/5 


4.4019256520 x 10" 


-7 


!J6 


3.9848995981 x 10" 


-8 


2/7 


5.3981503775 x 10" 


-15 



Table 9: Difference between starting 2/4 and constrained value 
m Difference 






-2.0767748211 x 10" 


-17 


1 


-2.0157837979 x 10" 


-17 


2 


-2.0157711737 x 10" 


-17 


3 


-2.0157697197 x 10" 


-17 


4 


-2.0157429383 x 10" 


-17 



from the starting value by no more than 10 and from each other by less, so are 
not particularly revealing to study directly. (Larger changes from the starting value 
could be obtained by starting from a different point with the same mass balance val- 
ues, but would not give any further insight.) Since it is difficult to compute the slow 
manifold (often a problem when real examples are used) we do not have a good way 
to characterize errors, but we can examine two features to see if the method appears 
to work. 

In Table ED we show the differences between the starting value and the constrained 
value of 2/4 for each order of constraint. We see that these differences show signs of 
"converging" - but this is certainly not irrefutable evidence of convergence. As a second 
test, we considered the relationship of the local derivative of the solution at the result 
of the constraint iteration. If it were very far from the slow manifold, we would expect 
it to have relatively large components of the eigenvectors corresponding to the large 
eigenvalues. (In general, even on the slow manifold it will not have zero components 
in the large eigendirections except for linear problems.) To estimate the amount of 
large eigencomponents present we computed v = dy/dt = f(y) and the local Jacobian 
J = df /dy at the solution, y, of the constraint iteration. Then we computed the norm 
ratio 




Its upper bound is the magnitude of the largest eigenvalue, and a value significantly 
less than this is an indicator that u is deficient in the largest eigendirection. Hence, 
the norm ratio gives some indication of the amount of the largest eigencomponents 
present. It values are shown in Table ITU1 

Since the largest eigenvalue is around 2.5 x 10 6 , it is clear that the m = case 
contains no more than around 20% of the large eigendirections, but this is drastically 
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Table 10: Norm Ratio ||Ji>||/||i>|| at Constrained Solution 

m R 

~0 4.44973316 x 10 5 

1 5.85785391 x 10 1 

2 6.18695075 x 10 1 

3 2.06227270 x 10 2 

4 1.50929245 x 10 2 



reduced for m = 1. From this particular starting point and choice of H, larger m gave 
no further improvement, but other choices of starting points or H yield norm ratios 
that reduce with each m increase, although by relatively small amounts. 



4.3 A Five-Dimensional System 

Because of the difficulty of determining whether the method is getting better ap- 
proximations as m increases, our final example is an artificial non-linear five-dimensional 
problem with a two-dimensional attractive invariant manifold. We start with the 
loosely coupled differential equations: 

dx\ 
~~dt ~ 

dx 2 



-X2 



Xl 



dw 2 2 \ 

— = L(x 1 +x 2 -w) 



~~dT 
du 2 
~~dt 



= (3 2 + u 2 2 



with L = 1000, f3i = 800, fo = 1200. The solutions of these are 

xi = Acos(t + cf)) 
X2 = Asin(t + (f)) 
w = A(l + be~ Lt ) 

Ui = -Pi/ii + ae-^) 

For any starting conditions, w — > x 2 (0) + x^O), and, if Uj(0) is chosen appropriately, 
Ui — > —Pi and the system goes to a closed orbit where the eigenvalues of the system 
Jacobian at each point on the closed orbit are =ti , -800, -100, and -1200. Thus w = 
x 2 + x\ is an attractive two dimensional invariant manifold. The above system is now 
subject to the unitary linear transformation given by y = Qv where v = [x T , w, u T ] T 
and Q is 



^ 5 



-3 
2 
2 
2 
2 



2 
-3 
2 
2 
2 
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Table 11: Residuals in Constraint Solutions. 



h 


m 


Residual 










Xo 




W 




U\ 








8.0 x icr 4 


1 


-4.84 x 10" 


-4 


-4.84 x 10" 


-4 


3.92 x 10" 


-3 


-2.50 x 10" 


-3 


-1.67 x 10" 


-3 




2 


-4.34 x 10" 


-6 


-4.34 x 10" 


-6 


2.55 x 10" 


-5 


-1.91 x 10- 


-5 


-8.50 x 10 


-6 




3 


-1.21 x 10" 


-6 


-1.21 x 10" 


-6 


-6.08 x 10" 


-7 


3.91 x 10- 


-9 


1.16 x 10 


-9 


2.0 x 10~ 4 


1 


-4.84 x 10" 


-4 


-4.84 x 10" 


-4 


3.92 x 10- 


-3 


-2.50 x 10" 


-3 


-1.67 x 10 


-3 




2 


-3.43 x 10" 


-6 


-3.43 x 10" 


-6 


2.59 x 10- 


-5 


-1.91 x 10- 


-5 


-8.50 x 10" 


-6 




3 


-3.01 x 10" 


-7 


-3.01 x 10" 


-7 


-1.55 x 10" 


-7 


3.73 x 10- 


-9 


1.14 x 10 


-9 


5.0 x 10~ 5 


1 


-4.84 x 10" 


-4 


-4.84 x 10" 


-4 


3.91 x 10- 


-3 


-2.49 x 10" 


-3 


-1.66 x 10 


-3 




2 


-3.22 x 10" 


-6 


-3.22 x 10" 


-6 


2.64 x 10- 


-5 


-1.95 x 10- 


-5 


-8.54 x 10" 


-6 




3 


-7.55 x 10" 


-8 


-7.55 x 10" 


-8 


-3.09 x 10" 


-8 


-7.27 x 10" 


-9 


4.20 x 10" 


10 



As in the first example, this is chosen to "mix up" the slow and fast components. We 
applied the constraint method using y\ and y2 as the fixed "observables." They were 
set to the values -791.2 and -792.2 respectively. The subspace y% = —791.2, 1/2 = —792.2 
intersects with the invariant manifold at four points (the defining system is a pair of 
quadratic equations). The intersection in the neighborhood of the solution has the 
values 

xi = -3.559434800714, x 2 = -2.559434800714 

Integration of y was performed using forward Euler with step size h for m + 1 steps, 
and iterating until the (m + l)-st forward difference was zero, with m = 1, 2, and 
3. The differences between the constraint calculations and the known solution are 
shown for several values of h in Table 1111 (They are called "residuals" there, since 
they are not exactly "errors.") We see that the residuals decrease by two to three 
orders of magnitude for each increase in m except for larger h. Larger h makes the 
iteration converge much more rapidly, but the increased error in the approximation 
of the difference to the derivative decreases the degree of approximation to the slow 
manifold. In a practical application, it might be wise to use a large h to get an initial 
approximation and then refine with a smaller h, although it might still be wise to use 
some convergence acceleration technique. 

5 Discussion and Conclusion 

We presented a "computational wrapper" approach for the approximation of a 
low-dimensional slow manifold using a legacy simulator. The approach effectively con- 
stitutes a protocol for the design and processing of brief computational experiments 
with the legacy simulator, which converge to an approximation of the slow manifold; in 
the spirit of CSP, one can think of it as "singular perturbation through computational 
experiments". It is interesting that, if one could initialize a laboratory experiment 
at will, our "computer experiment" protocol could become a laboratory experiment 
protocol for the approximation of a slow manifold. 
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The approach can be enhanced in many ways; we already mentioned the possible 
use of matrix-free fixed point algorithms for the acceleration of its convergence. Here 
we used the "simplest possible" estimation (through finite difference interpolation) of 
the trajectory from the results of the simulation. Better estimation techniques (e.g. 
maximum likelihood) can be linked with the data processing part of the approach; this 
will be particularly important when the results of the detailed integration are noisy, 
as will be the case in the observation of the evolution of statistics of complex evolving 
distributions. 

It is also important to notice that, upon convergence of the procedure, one can 
implement a matrix-free, timestepper based computational approximation of the lead- 
ing eigenvalues of the local linearization of the dynamics (e.g. through a timestepper 
based Arnoldi procedure, see [5,29]). As the evolution progresses, or as the parameters 
change, this test can be used to adaptively adjust the local dimension of the slow man- 
ifold - we can detect whether a slow mode is starting to become fast, or when a mode 
that used to be fast is now becoming slow. The eigenvectors of these modes constitute 
good additional observables for the parameterization of a "fatter" slow manifold. One of 
the important features of the approach is that one does not need to a priori know what 
the so-called "slow variables" are - any set of observables that can parameterize the 
slow manifold (i.e., over which the manifold is the graph of a function) can be used for 
our approach. If data analysis ( [14,30]) suggests good observables that are nonlinear 
combinations of the "obvious" state variables, the approach can still be implemented; 
the knowledge of good "order parameters" can thus be naturally incorporated in this 
approach. 

Overall, this approach provides us with a good initial condition for the full problem, 
consistent with a set of observables - an initial condition that lies close to the slow 
manifold, sometimes referred to as a "mature" or "bred" initial condition. Such initial 
conditions are essential for the implementation of equation- free algorithms: algorithms 
that solve the reduced problem without ever deriving it in closed form [4, 25, 27]). 
Indeed, short bursts of appropriately initialized simulations can be used to perform long 
term prediction (projective and coarse projective integration) for the reduced problem, 
its stability and bifurcation analysis, as well as tasks like control and optimization. 
We expect this approach to become a vital component of the "lifting" operator in 
equation-free computation. 
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